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The  response  of  the  energetic  molecular  crystal  cyclotrimethylene  trinitramine  (RDX)  to  the  propagation  of 
planar  shock  waves  normal  to  (100)  has  been  studied  using  large-scale  molecular  dynamics  simulations  that 
employ  an  accurate  and  transferable  nonreactive  potential.  The  propagation  of  the  shock  waves  was  simulated 
using  nonequilibrium  molecular  dynamics.  Shear  bands  were  nucleated  during  shocks  with  a  particle  velocity 
of  1.0  km  s_1  and  corresponding  Rankine-Hugoniot  shock  pressure  of  9.7  GPa.  These  defects  propagate  into 
the  compressed  material  at  45°  to  [100]  in  the  [010]  zone.  The  shear  bands  evolve  slowly  compared  to  the  time 
scales  routinely  accessible  to  nonequilibrium  molecular  dynamics  toward  a  liquidlike  state  as  a  result  of 
viscous  heating.  A  recently  developed  shock-front  absorbing  boundary  condition  [A.  V.  Bolesta  et  al. ,  Phys. 

Rev.  B  76,  224108  (2007)]  was  applied  to  the  simulation  cells  at  the  moment  of  maximum  compression  to 
sustain  the  shock-compressed  state.  Molecular  dynamics  simulations  were  then  employed  to  study  the  temporal 
and  structural  evolution  of  the  shock-induced  shear  bands  toward  a  steady-fluctuating  state.  Owing  to  the 
intense,  viscous  flow-driven  heating  within  the  shear  bands,  these  defects  can  be  considered  to  be  homoge¬ 
neously  nucleated  hot  spots. 

DOI:  10.1103/PhysRevB.78.014107  PACS  number(s):  62.50.-p,  62.20.F-,  61.72.Bb,  61.20.Ja 


I.  INTRODUCTION 

Defects  in  solid  and  crystalline  energetic  materials  are 
known  to  exert  a  significant  influence  on  impact  and  initia¬ 
tion  sensitivity.1  Mesoscopic  and  macroscopic  defects  lead  to 
the  spatial  localization  of  the  translational  energy  from  a 
shock  through,  for  example,  the  formation  of  jets  during  the 
collapse  of  voids2,3  or  interfacial  friction  at  cracks  or  particle 
boundaries.4,5  Events  of  this  kind  result  in  hot  spots  where 
temperature  and/or  stress  may  exceed  significantly  values  in 
the  bulk  and  hence  promote  molecular  decomposition.  The 
formation  of  hot  spots  through  these  mechanisms  is  well 
understood,  at  least  on  a  qualitative  level,6  and  apply  in  gen¬ 
eral  to  all  solid  energetic  materials. 

The  impact  sensitivity  of  defect-free  single  crystals  of  the 
energetic  molecular  crystal  pentaerythritol  tetranitrate 
(PETN)  was  found  to  depend  strongly  on  the  crystallo¬ 
graphic  orientation  of  the  shock  propagation  direction.  7~y  In 
this  case,  dislocation-mediated  plastic  deformation  was  pro¬ 
posed  as  an  explanation  for  the  observed  orientation  depen¬ 
dencies.  Slip  systems  were  identified  in  single  crystals  using 
x-ray  topography  and  by  the  analysis  of  slip  traces  after  sur¬ 
face  indentation.10  Those  orientations  of  the  shock  propaga¬ 
tion  direction  for  which  there  was  no  resolved  shear  stress  on 
any  slip  system  were  found  to  have  a  high  Hugoniot  elastic 
limit  and  high  impact  sensitivity.  Similarly,  those  orienta¬ 
tions  for  which  a  shear  stress  was  resolved  on  all  of  the  slip 
systems  identified  experimentally  were  found  to  have  low 
Hugoniot  elastic  limit  and  low  impact  sensitivity.  Hence,  in 
the  absence  of  mesoscopic  and  macroscopic  defects,  ex¬ 
tended  crystal  defects  at  the  atomic  or  molecular  scale  were 
found  to  play  a  role  in  determining  impact  sensitivity.  This 
work  led  to  the  formulation  of  the  steric-hindrance 


model,7,8,11  wherein  an  absence  of  dislocation-mediated  plas¬ 
ticity  for  a  given  crystallographic  orientation  of  the  shock 
propagation  direction  promotes  initiation  since  shear  stresses 
cannot  be  relaxed  easily  and  molecules  deform  severely,  in¬ 
ducing  bond-breaking  events. 

We  have  studied  the  response  of  single  crystals  of  the 
widely  used  energetic  molecular  crystal  cyclotrimethylene 
trinitramine,  C3H6N606  (RDX),  to  the  propagation  of  shock 
waves  normal  to  (100)  by  means  of  large-scale  molecular 
dynamics  (MD)  simulations,  with  particular  focus  on  the 
shock-induced  nucleation  of  extended  defects.  Under  ambi¬ 
ent  pressure  and  temperature,  RDX  adopts  an  orthorhombic 
unit  cell  in  space  group  Pbca  that  contains  eight  molecules 
(o'-polymorph).12  An  RDX  molecule  is  depicted  in  Fig.  1(a) 
and  a  projection  of  the  a-RDX  unit  cell  along  [001]  in  Fig. 
1(b).  Three  slip  systems  were  identified  in  RDX  single  crys¬ 
tals  by  the  analysis  of  slip  traces  after  surface  indentation: 
(010)[001],  (02 1)[100],  and  (02l)[100].10  Hence,  the  [100] 
shock  propagation  direction  is  special  since  there  is  no  re¬ 
solved  shear  stress  on  any  of  the  slip  systems  identified  un- 


FIG.  1.  (Color  online)  (a)  An  RDX  molecule  and  (b)  projection 
of  the  a-RDX  unit  cell  along  [001]. 
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der  quasistatic  loading.13  Thus,  on  the  basis  of  the  steric- 
hindrance  model,  we  expect  this  orientation  to  be  of  high 
impact  sensitivity.  However,  earlier  work  on  (lll)-oriented 
RDX  single  crystals  revealed  an  abrupt  change  in  the  mecha¬ 
nism  of  plastic  deformation  above  a  certain  shock  pressure 
that  could  not  have  been  anticipated  based  on  either  crystal¬ 
lographic  arguments  or  the  results  of  quasistatic  loading 
experiments.1314  This  work  led  us  to  the  conclusion  that  ex¬ 
trapolating  the  mechanisms  for  plastic  deformation  identified 
in  energetic  molecular  crystals  under  quasistatic  loading  over 
many  orders  of  magnitude  in  strain  rate  can  be  rather  unre¬ 
liable,  particularly  for  complicated,  low-symmetry  crystal 
structures.  In  fact,  PETN  can  be  considered  as  a  special  case 
in  this  regard  since  it  adopts  a  relatively  high-symmetry  te¬ 
tragonal  unit  cell  comprising  only  two  molecules. 

Nonequilibrium  molecular  dynamics  (NEMD) 
simulations15"20  of  the  propagation  of  planar  shock  waves 
normal  to  (100)  were  performed  for  particle  velocities  U p 
=  0.63  and  1.0  km  s_1,  corresponding  to  shock  pressures, 
PRH,  calculated  using  the  Rankine-Hugoniot  jump 
conditions21  of  5.7  and  9.7  GPa,  respectively.  At  PRH 
=  5.7  GPa,  we  found  no  evidence  for  the  nucleation  of  crys¬ 
tal  defects.  However,  at  PRH=9.7  GPa,  liquidlike  shear 
bands  (SBs)  were  nucleated  which  propagate  at  45°  to  the 
compression  axis.  The  remainder  of  this  paper  is  devoted  to 
the  characterization  of  these  defects. 

A  popular  method  for  generating  planar  shock  waves  in 
NEMD  simulations  is  to  either  drive  a  rigid  piston  at  a  speci¬ 
fied  particle  velocity  onto  a  stationary  simulation  cell  or  im¬ 
pact  a  simulation  cell  onto  a  fixed,  rigid  piston  at  a  specified 
particle  velocity.  In  this  paper  we  have  employed  the  latter 
approach.  Both  of  these  methods  limit  the  time  interval  over 
which  material  remains  on  the  Hugoniot  locus  since  once  the 
shock  wave  reaches  the  free  surface  of  the  simulation  cell,  a 
rarefaction  wave  propagates  rapidly  back  into  the  shock- 
compressed  material.  Hence,  material  near  the  free  surface  is 
on  the  Hugoniot  locus  for  a  very  limited  time.  Simulations  of 
this  type  typically  employ  simulation  cells  that  are  relatively 
long  parallel  to  the  shock  propagation  direction  in  order  to 
maximize  the  time  interval  over  which  material  in  the  vicin¬ 
ity  of  the  piston  is  sustained  in  the  shocked  state.  This  ap¬ 
proach  is  often  not  practical  computationally,  particularly  in 
the  simulation  of  slow  processes  such  as  plasticity  or  chem¬ 
istry.  In  the  present  work,  this  limitation  of  NEMD  simula¬ 
tion  techniques  is  particularly  pronounced  since  the  internal 
structure  of  the  shock-induced  shear  bands  evolves  relatively 
slowly,  compared  to  typical  NEMD  time  scales,  by  a  viscous 
heating  mechanism. 

We  have  employed,  with  modifications,  a  simple  and  ro¬ 
bust  method  for  extending  time  scales  in  molecular  dynamics 
simulations  of  shock  loading  that  was  first  developed  by 
Bolesta  et  al.22  The  shock-front  absorbing  boundary  condi¬ 
tion  (SFABC)  (Ref.  23)  enables  a  seamless  transition  from  a 
NEMD  simulation  of  the  propagation  of  a  shock  front  to  the 
simulation  of  shock-induced  defects  during  their  evolution 
toward  a  steady-fluctuating  state.  SFABCs  capture  the  simu¬ 
lation  cell  at  the  point  of  maximum  compression  in  a  NEMD 
simulation  and  prevent  the  emission  of  rarefaction  waves 
from  the  free  surface.  Furthermore,  the  SFABC  approach  ob¬ 
viates  the  requirement  for  simulation  cells  that  are  “long” 


parallel  to  the  direction  of  shock  propagation,  does  not  intro¬ 
duce  any  incoherent  interfaces  into  the  system,  and  leaves 
unaffected  the  microcanonical  equations  of  motion.  The 
SFABC  method  and  its  application  to  a  crystalline  solid  are 
described  in  detail  in  Sec.  III. 

Other  methods  for  the  absorption  of  waves  incident  at 
boundaries  in  MD  simulations  have  been  developed  in  recent 
years.  For  example,  Namilae  et  al.24  implemented  a  differen¬ 
tial  equation-based  absorbing  boundary  condition  (ABC)  to 
match  the  impedance  of  a  semi-infinite  continuum  space  to 
that  of  an  atomistic  region.  However,  ABCs  of  this  type  are 
based  on  the  linear  wave  equation  and  for  this  reason  are  not 
suited  to  the  absorption  of  shock  waves. 

H.  NONEQUILIBRIUM  MOLECULAR  DYNAMICS 
SIMULATIONS 

The  nonreactive,  fully  flexible  molecular  potential  for  nit- 
ramines  developed  by  Smith  and  co-workers25-26  was  em¬ 
ployed  in  all  of  the  MD  simulations.  The  Smith-Bharadwaj 
potential25  describes  bond  stretches,  bond  angles,  and  out-of- 
plane  bends  using  harmonic  springs.  Torsions  are  represented 
by  anharmonic  terms  that  display  extrema  at  the  torsion 
angles  that  correspond  to  stationary  points  on  the  conforma¬ 
tional  energy  surface.  The  intramolecular  terms  were  param¬ 
etrized  to  quantum  chemistry  calculations  of  the  structure, 
vibrational  frequencies,  and  barriers  to  conformational 
change  in  a  model  compound,  1, 3-dimethyl- 1,3-dinitro  me- 
thyldiamine.  Intermolecular  interactions  are  represented  by  a 
sum  of  Buckingham  potentials  parametrized  to  standard  lit¬ 
erature  values.  Electrostatic  interactions  are  included  explic¬ 
itly  and  employ  partial  charges  that  were  increased  for  all 
species  by  25%  compared  to  quantum  chemistry-calculated 
values  to  account  for  the  effects  of  polarization  in  condensed 
phases.  The  Smith-Bharadwaj  potential,25  while  not  param¬ 
etrized  to  reproduce  any  particular  property  of  condensed 
phase  nitramines,  was  shown  to  provide  excellent  descrip¬ 
tions  of  crystalline  cyclotetramethylene  tetranitramine 
(HMX)  (Refs.  26  and  27)  and  RDX.  In  the  case  of  RDX  the 
Smith-Bharadwaj  potential25  not  only  predicts  the  ortho¬ 
rhombic  Pbca  space  group  to  be  the  crystal  structure  with 
the  lowest  free  energy  at  zero  pressure  and  300  K  but  it  also 
provides  remarkably  accurate  predictions  for  its  lattice 
parameters,12  single-crystal  elastic  constants,28  coefficient  of 
thermal  expansion,29  and  isothermal  compression  curve  up  to 
the  a-to-y  phase  transformation.30  Furthermore,  it  predicts 
the  P  versus  U p  Hugoniot  of  RDX  crystals  shock  loaded 
normal  to  (111)  to  within  0.5  GPa  of  experiment.31 

The  lattice  parameters  of  a-RDX  calculated  using  the 
Smith-Bharadwaj  potential25  at  300  K  and  zero  pressure  are 
a  =  13.400  A,  b=  11.517  A,  and  c=  10.623  A,  each  of 
which  is  within  1.7%  of  experiment. 12  These  calculated  lat¬ 
tice  parameters  were  used  for  all  the  simulations  reported 
here.  NEMD  simulations  of  the  propagation  of  planar  shock 
waves  normal  to  (100)  employed  an  oblong  simulation  cell 
measuring  120a  X  20b  X  20c,  containing  8.064X  106  atoms. 
A  shorter  simulation  cell  measuring  60a  X  20 b  X  20c  was 
also  used  to  study  the  effect  of  SFABCs  on  the  structure  of 
shock-induced  defects.  These  cells  are  hereafter  referred  to 
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FIG.  2.  (a)  Schematic  illustration  of  a  NEMD  shock  loading 
simulation.  The  vertical  broken  line  represents  the  position  of  the 
shock  front  which  is  traveling  from  left  to  right,  (b)  Simulation  cell 
at  maximum  compression,  (c)  Application  of  the  second  piston  in 
the  shock-front  absorbing  boundary  condition  (SFABC)  method. 

as  the  long  and  short  cells,  respectively.  Finally,  we  also 
created  a  quasi-two-dimensional  (quasi-2D)  simulation  cell 
to  assess  and  mitigate  any  finite-size  effects  on  the  defect 
structures  generated  during  the  simulations.  The  quasi-2D 
cell  measured  149a  X  3 bX  141  c  and  contained  10.6  X  106  at¬ 
oms.  An  additional  200  A  thick  vacuum  region  was  included 
in  the  simulation  cells  along  [100]  to  minimize  electrostatic 
interactions  between  free  surfaces  upon  the  application  of 
three-dimensional  periodic  boundary  conditions.  The  simula¬ 
tion  cells  were  thermalized  until  intermolecular  and  intramo¬ 
lecular  temperatures  equilibrated  to  300  K  prior  to  shock 
loading.  The  propagation  of  supported  planar  shock  waves 
was  simulated  in  the  microcanonical  ensemble  by  impacting 
the  cells  onto  a  fixed  piston  of  thickness  3a  consisting  of 
rigid  RDX  molecules  in  the  same  crystallographic  orienta¬ 
tion  as  the  simulation  cell  by  adding  to  all  atoms  a  particle 
velocity  U p  parallel  to  [100].  This  generates  a  shock  wave 
propagating  at  velocity  Uw=Us—Up  relative  to  the  stationary 
piston,  where  Us  is  the  shock  wave  velocity;  this  is  illus¬ 
trated  schematically  in  Fig.  2(a).  All  simulations  were  per¬ 
formed  using  the  LAMMPS  code.32  Long-range  electrostatic 
interactions  were  calculated  using  the  PPPM  method33  and  all 
C-H  bonds  were  constrained  to  equilibrium  length  using  the 
SHAKE  algorithm.  One  NEMD  simulation  was  performed  at 
Up=0.63  km  s-1  using  the  long  cell  and  three  at  Up 
=  1.0  km  s-1  using  the  long,  short,  and  quasi-2D  cells.  A  0.4 
fs  time  step  for  the  integration  of  the  equations  of  motion 
was  employed  in  each  case. 

III.  SHOCK-FRONT  ABSORBING  BOUNDARY 
CONDITIONS 

The  SFABC  method  facilitates  a  seamless  transition  from 
a  NEMD  simulation  of  the  propagation  of  a  shock  wave 
through  a  medium  to  a  MD  simulation  of  the  evolution  of 
any  defects  created  by  the  shock  wave.  SFABCs  were  origi¬ 
nally  implemented  and  applied  to  MD  studies  of  shock 
waves  in  reactive  methane.22  In  the  present  work  on  crystal¬ 
line  RDX,  we  used  a  modest  simplification  of  the  original 
method  for  absorbing  the  shock  front  which  we  expect  to  be 
more  robust  in  its  implementation.  The  SFABC  method  is 
illustrated  schematically  in  Fig.  2. 
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FIG.  3.  Plots  of  the  x  coordinate  of  the  last  molecule  as  a  func¬ 
tion  of  the  time  step  from  the  NEMD  simulations  of  (a)  the  long 
cell,  Up  =  0.63  km  s-1,  (b)  long  cell,  Up=  1.0  km  s-1,  and  (c)  short 
cell,  Up  =1.0  kms-1. 

The  starting  point  for  the  application  of  SFABCs  is  a 
NEMD  simulation  of  shock  loading  as  described  in  Sec.  II. 
The  shock  wave  is  allowed  to  propagate  until  the  simulation 
cell  reaches  the  point  of  maximum  compression.  The  point  of 
maximum  compression  can  be  identified  to  a  high  level  of 
accuracy  by  monitoring  as  a  function  of  time  the  x  coordi¬ 
nate  of  the  center  of  mass  of  the  “last”  molecule  in  the  simu¬ 
lation  cell.  The  last  molecule,  the  one  with  the  largest  x  co¬ 
ordinate,  can  be  identified  easily  by  postprocessing  analysis 
of  the  NEMD  simulation.  The  simulation  cells  were  captured 
at  maximum  compression  to  within  a  tolerance  of  1  A  in  the 
three  NEMD  simulations  described  in  Sec.  II.  Plots  of  the  x 
coordinate  of  the  last  molecule  as  a  function  of  the  time  step 
are  presented  in  Fig.  3. 

Once  the  point  of  maximum  compression  was  identified 
in  the  NEMD  simulations,  the  force  and  velocity  vectors  of 
all  atoms  within  the  last  30  A  of  the  cell  were  set  equal  to 
zero.  In  this  manner,  those  molecules  at  the  end  of  the  simu¬ 
lation  cell  are  treated  as  a  second  piston  that,  together  with 
the  first  stationary  rigid  piston,  serve  to  constrain  the  shock- 
compressed  material  at  constant  volume.  In  principle,  it 
would  also  be  possible  to  capture  the  simulation  cell  at  maxi¬ 
mum  compression  by  removing  from  the  system  the  vacuum 
region  and  applying  periodic  boundary  conditions  along  the 
shock  propagation  direction.  However,  the  approach  adopted 
here  ensures  that  the  interface  between  the  simulation  cell 
and  second  piston  is  structurally  coherent,  unlike  the  inter¬ 
face  that  would  be  created  upon  the  application  of  periodic 
boundary  conditions.  Furthermore,  since  the  second  piston 
can  be  made  arbitrarily  thick,  the  SFABC  method  is  not  af¬ 
fected  adversely  by  the  finite  width  of  the  shock  front. 

In  the  original  implementation  of  the  SFABC  method,22 
the  point  of  maximum  compression  was  determined  by 
monitoring  as  a  function  of  time  the  specific  energy  content 
of  layers  of  molecules  perpendicular  to  the  direction  of  shock 
propagation.  The  profile  of  specific  energy  versus  position 
was  fit  to  a  straight  line  which  was  extrapolated  to  predict  the 
time  at  which  the  second  piston,  which  is  this  case  was  rigid 
and  at  rest,  was  assigned  velocity  Up. 
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FIG.  4.  Projections  along  [010]  of  the  molecular  centers  of  mass 
after  shock  loading  at  Up=  1.0  km  s-1.  (a)  Long  cell  at  maximum 
compression  (Ref.  34),  (b)  long  cell  after  121  ps  MD  simulation 
under  SFABCs  (Ref.  34),  (c)  short  cell  at  maximum  compression, 
and  (d)  short  cell  after  110  ps  MD  simulation  under  SFABCs.  The 
crystalline  regions  (CRs)  and  shear  bands  (SBs)  used  for  subse¬ 
quent  analysis  are  indicated  in  (b)  and  (d). 

IV.  RESULTS 
A.  Shock-induced  defects 

1.  NEMD  simulations 

The  NEMD  simulation  of  the  propagation  of  a  planar 
shock  wave  along  [100]  in  the  long  cell  at  Up 
=  0.63  kms"1  (Prh=5.7  GPa)  did  not  result  in  the  nucle- 
ation  of  crystal  defects.  As  noted  in  Sec.  I,  at  this  orientation 
of  the  shock  propagation  direction,  none  of  the  dislocation 
slip  systems  identified  under  quasistatic  loading  conditions10 
are  subject  to  a  resolved  shear  stress. 

Shock  loading  at  Up=l.O  kms-1  (PRH=9.7  GPa)  led  to 
the  formation  of  shear  bands  during  NEMD  simulations  em¬ 
ploying  both  the  long  and  short  cells.  Projections  of  the  mo¬ 
lecular  centers  of  mass  along  [010]  of  the  long  and  short 
cells  at  maximum  compression  are  presented  in  Figs.  4(a) 
and  4(c),  respectively.35  In  both  cases,  the  shock  wave  propa¬ 
gated  from  left  to  right  and  it  is  apparent  that  the  shear  bands 
are  structurally  more  developed  in  the  vicinity  of  the  piston 
since  their  growth  rate  is  relatively  slow  compared  to  the 
time  scale  of  the  simulations.  The  shear  bands  propagate  at 
approximately  45°  to  [100]  and  the  plane  of  the  bands  lies  in 
the  [010]  zone.  There  appears  to  be  no  preference  as  to 
whether  the  shear  bands  propagate  at  positive  or  negative 
45°  to  [100];  the  bands  are  not  aligned  with  any  crystallo¬ 
graphic  orientation  of  the  underlying  structure  and  are  driven 
only  by  shear  stresses  which  are  a  maximum  at  these  orien¬ 
tations.  However,  there  is  a  distinct  preference  for  the  plane 
of  the  bands  to  lie  in  the  [010]  zone.  The  shear  bands  are  not 
appreciably  higher  in  density  than  the  surrounding  crystalline 
regions  (CRs);  the  latter  appear  to  be  of  low  density  since  in 
this  projection  columns  of  molecules  are  viewed  “end  on,” 
while  the  former  comprise  amorphous  material. 

2.  SFABC  simulations 

The  SFABC  was  applied  to  the  simulation  cells  as  de¬ 
scribed  in  Sec.  III.  The  long  cell  shock  loaded  at  Up 


=  0.63  km  s-1  was  simulated  in  the  microcanonical  ensemble 
for  42  ps  after  the  application  of  SFABCs.  The  long  and 
short  cells  shock  loaded  at  Up=  1.0  km  s_1  were  simulated  in 
the  microcanonical  ensemble  under  SFABCs  for  121  and  110 
ps,  respectively.  Snapshots  of  the  long  and  short  cells  shock 
loaded  at  Up=  1 .0  kms-1  taken  at  the  conclusion  of  the 
SFABC  simulations  are  presented  in  Figs.  4(b)  and  4(d),  re¬ 
spectively.  The  SFABCs  enable  simulation  times  sufficient 
for  the  shear  bands  to  propagate  throughout  the  entire  sys¬ 
tem,  and,  more  importantly,  for  their  internal  structure  to 
evolve  to  a  steady-fluctuating  state.  Detailed  analyses  of  the 
internal  structure  of  the  shear  bands  are  provided  in  Sec. 
IV  B. 

Upon  the  application  of  SFABCs,  shear  bands  propagate 
into  the  shock-compressed  material  from  the  second  piston 
in  the  opposite  orientation  to  those  formed  at  the  first  piston. 
From  the  moment  that  the  second  piston  is  applied  and  both 
ends  of  the  cell  are  fixed,  the  total  plastic  strain  mediated  by 
the  shear  bands  in  the  directions  normal  to  [100]  must  equal 
to  zero.  Hence,  growth  of  the  shear  bands  nucleated  in  the 
NEMD  simulation  is  compensated  by  shear  in  the  opposite 
sense  provided  by  shear  bands  growing  from  the  opposite 
end  of  the  cell.  While  this  may  at  first  seem  to  be  an  unphysi¬ 
cal  artifact  associated  with  SFABCs,  it  is,  in  fact,  physically 
correct.  We  previously  observed  an  abrupt  change  in  the  ori¬ 
entation  of  shear  bands  in  NEMD  simulations  of  systems 
with  very  high  aspect  ratio  (15:1  rather  than  8:1)  owing  to 
the  constraints  on  lateral  motion  imposed  by  the  unshocked 
material  ahead  of  the  shock  front.36  Hence,  in  this  respect, 
SFABCs  mimic  a  NEMD  simulation  using  a  cell  of  very  high 
aspect  ratio. 

Simulations  using  the  long  and  short  cells  at  /JRH 
=  9.7  GPa  (Up=l  km  s-1)  exhibit  pronounced  finite-size  ef¬ 
fects  whereby  one  shear  band  propagates  a  significant  dis¬ 
tance  through  the  simulation  cell  owing  to  application  of 
periodic  boundary  conditions.  While  such  finite-size  effects 
are  generally  undesirable,  they  in  no  way  affect  our  conclu¬ 
sions.  In  Fig.  5  we  present  a  snapshot  of  the  quasi-2D  cell 
after  shock  loading  to  9.7  GPa  using  NEMD  and  56.8  ps  of 
MD  simulation  after  the  application  of  the  SFABC.  The  ori¬ 
entation  and  dimensions  of  the  quasi-2D  cell  were  selected 
both  to  lead  to  the  formation  of  shear  bands  (the  shortest  axis 
of  the  cell  is  parallel  to  [010])  and  to  minimize  the  role  of 
finite-size  effects  on  their  subsequent  growth.  It  is  clear  from 
Fig.  5  that  many  shear  bands  have  nucleated  randomly  both 
at  the  pistons  and  in  the  bulk.  Furthermore,  the  shear  bands 
have  grown  at  ±45°  to  [100]  in  roughly  equal  amounts. 
Since  only  a  small  fraction  of  the  shear  bands  cross  the  (001) 
periodic  boundaries,  we  propose  that  the  simulated  spatial 
distribution  of  defects  corresponds  closely  to  that  which 
would  be  observed  in  a  mesoscopic  specimen. 

B.  Structural  and  thermodynamic  analysis 

The  simulations  at  Up=0.63  and  1.0  km  s_l  were  exam¬ 
ined  to  quantify  shock-induced  changes  in  molecular  confor¬ 
mation,  the  internal  structure  of  the  shear  bands,  and  the 
temperature  rise  caused  by  viscous  heating  during  their 
nucleation  and  growth.  We  studied  these  quantities  both  for 
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FIG.  5.  A  projection  of  the  molecular  centers  of  mass  of  the 
quasi-2D  simulation  cell  following  shock  loading  to  9.7  GPa  using 
NEMD  and  a  56.8  ps  MD  simulation  after  the  application  of  the 
SFABC  (Ref.  34).  At  maximum  compression,  the  cell  measures 
1638.5X34.55X  1497.9  A3  parallel  to  [100],  [010],  and  [001], 
respectively. 

crystalline  regions  and  shear  bands  adjacent  to  these  regions 
in  our  simulations  at  Up=l.O  km  s-1.  The  sampled  CRs  and 
SBs  are  indicated  in  Figs.  4(b)  and  4(d). 

1.  Radial  distribution  functions 

Radial  distribution  functions  (RDFs)  for  molecular  cen¬ 
ters  of  mass  were  calculated  to  determine  whether  the  shear 
bands  comprise  a  liquidlike  structure  that  is  consistent  with 
melting  and/or  amorphization.  The  RDF  is  defined  as  g(R) 
=  n(R)l4TrpR2AR ,  where  ii(R)  is  the  number  of  particles  in  a 
spherical  shell  of  radius  R  and  thickness  A R  and  p  the  par¬ 
ticle  number  density.  We  used  A// =  0.2  A  in  all  calculations 
and  p  was  determined  uniquely  for  each  system.  All  RDFs 
were  averaged  over  at  least  7900  molecules. 

We  plot  the  RDFs  for  an  unshocked  RDX  crystal  at  300  K 
and  a  crystalline  region  at  the  point  of  maximum  compres¬ 
sion  from  a  NEMD  simulation  at  Up=  1.0  km  s'1  in  Figs. 
6(a)  and  6(b),  respectively.  Comparing  Figs.  6(a)  and  6(b), 
we  find  only  a  small  change  in  the  RDFs  at  first  nearest 
neighbors  but  notable  changes  at  second-nearest  neighbors 
and  beyond.  The  most  relevant  feature  is  the  decreased  gap 
between  first-  and  second-nearest  neighbors.  This  is  expected 
due  to  the  large  uniaxial  compression  imparted  by  shock 
loading.  The  RDF  of  the  shock-compressed  system  becomes 
relatively  smooth  at  R>16  A  which  suggests  some  loss  in 
long-range  order,  although  the  system  is  still  clearly  crystal¬ 
line. 

The  RDFs  calculated  for  a  crystalline  region  and  an  adja¬ 
cent  shear  band  after  a  121  ps  MD  simulation  under  SFABCs 
are  presented  in  Figs.  7(a)  and  7(b),  respectively.  The  RDF 
for  the  crystalline  region  shows  a  recovery  of  the  gap  be¬ 
tween  the  first-  and  second-nearest-neighbor  shells  that  was 


FIG.  6.  Center-of-mass  radial  distribution  functions,  (a)  Un¬ 
shocked  a-RDX  crystal  at  300  K  and  (b)  crystalline  region  from  the 
long  simulation  cell  at  maximum  compression  after  shock  loading 
at  Up=  1.0  km  s-1. 

reduced  during  the  NEMD  phase  of  the  simulation  [Fig. 
6(b)].  This  suggests  that  SFABCs  facilitate  the  evolution  of 
the  system  from  a  metastable  configuration  after  shock  load¬ 
ing  to  geometries  more  similar  to  those  in  the  unshocked 
crystal.  The  RDF  for  the  shear  band  corresponds  unambigu¬ 
ously  to  that  of  a  liquidlike  state  and  shows  that  the  shear 
bands  are  regions  of  localized,  shear  stress-driven  amor¬ 
phization. 

2.  Heating  via  viscous  flow 

The  intermolecular  and  intramolecular  temperatures  of  a 
crystalline  region  and  adjacent  shear  band  were  calculated  as 
a  function  of  time  during  equilibrium  MD  simulations  em¬ 
ploying  SFABCs.  Plots  of  the  temperatures  calculated  in  the 
long  and  short  cells  are  shown  in  Figs.  8  and  9,  respectively. 
In  both  the  long  and  short  cells,  the  onset  of  shear  band 
formation  can  be  discerned  clearly  by  the  rapid  increase  in 
the  intermolecular  and  intramolecular  temperatures  that  is 


FIG.  7.  Center-of-mass  radial  distribution  functions  after  121  ps 
MD  simulation  under  SFABCs.  (a)  Crystalline  region  and  (b)  shear 
band. 
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FIG.  8.  Intermolecular  and  intramolecular  temperatures  of  a 
crystalline  region  and  adjacent  shear  band  during  an  MD  simulation 
employing  SFABCs  and  the  long  cell.  Zero  time  corresponds  to  the 
start  of  the  SFABC  simulation. 

driven  by  viscous  flow.34  In  both  simulations,  the  intramo¬ 
lecular  temperature  within  the  shear  band  increases  from 
around  390  to  a  maximum  of  460  K.  The  passage  of  the 
shock  wave  is  responsible  for  the  90  K  increase  in  intramo¬ 
lecular  temperature  compared  to  that  in  the  unshocked  RDX 
crystal.  The  heating  rate  within  the  shear  band  is  about  3.5 
X  1012  K  s-1  during  the  initial  stage  of  growth  in  the  long 
cell  and  about  5.8  X  1012  K  s-1  in  the  short  cell. 

The  intense  direct  heating  of  material  in  the  shear  bands 
during  viscous  flow  drives  the  indirect  heating  of  the  sur¬ 
rounding  crystalline  regions.  After  the  initial  growth  of  the 
shear  bands,  the  intermolecular  and  intramolecular  tempera¬ 
tures  equilibrate  throughout  the  crystal.  Both  simulations 
suggest  that  temperatures  of  the  shear  bands  and  bulk  will 
equilibrate  to  about  430  K.  Here  we  are  using  a  nonreactive 
potential,  but  in  reality  the  intense  heating  and  viscous  flow 
associated  with  the  formation  of  shear  bands  would  cause 
some  RDX  molecules  to  decompose  via  a  strongly  exother¬ 


mic  path.  Assuming  that  this  chemistry  progresses  to  the 
exothermic  steps  in  the  decomposition  mechanism,  the  heat¬ 
ing  would  be  enhanced  by  the  processes.  Nevertheless,  the 
release  of  stored  elastic  strain  energy  that  the  formation  of 
shear  bands  facilitates  increases  the  temperature  of  the  entire 
crystal  by  40-50  K  over  that  due  to  the  passage  of  the  shock 
wave  alone. 

It  is  important  to  note  that  the  levels  of  shock  heating  we 
calculate  from  our  classical  MD  simulations  underestimate 
those  that  would  take  place  in  a  real  organic  molecular  crys¬ 
tal.  Within  classical  MD,  each  degree  of  vibrational  freedom 
is  associated  with  a  thermal  energy  kBTI  2.  However,  the  De¬ 
bye  temperatures  of  the  intramolecular  vibrational  modes  in 
RDX  are  known  to  be  significantly  higher  than  room  tem¬ 
perature  (see,  for  example.  Ref.  37).  Thus,  the  true  heat  ca¬ 
pacity  of  RDX  under  ambient  conditions  is  significantly  less 
than  the  classical  limit. 

The  shock  compression  determined  from  our  NEMD 
simulations  of  shock  loading  to  PRH=9.7  GPa  is  V/V0 
=  0.82.  This  can  be  considered  a  weak  shock  since  PRH  is  less 
than  the  bulk  modulus  of  the  material.  Using  experimental 
and  experimentally  derived  values  for  material  properties, 
we  can  estimate  the  shock  heating  that  takes  place  in  real 
RDX  under  the  same  loading  conditions  by  summing  contri¬ 
butions  from  the  heating  upon  isentropic  compression  from 
specific  volume  V0  to  V,  A Ts,  and  from  the  additional  heat¬ 
ing  that  arises  from  the  passage  of  the  shock  front,  ATH, 


at=ats+ath. 


(1) 


where 


A  TS=T0 


(2) 


and 


FIG.  9.  Intermolecular  and  intramolecular  temperatures  of  a 
crystalline  region  and  adjacent  shear  band  during  an  MD  simulation 
employing  SFABCs  and  the  short  cell.  Zero  time  corresponds  to  the 
start  of  the  SFABC  simulation. 


AT, 


"J. 


de 

es(V)  Cv 


(3) 


Here,  T  is  the  Griineisen  coefficient,  T0  the  initial  tempera¬ 
ture,  cv  the  heat  capacity  at  constant  volume,  es(  V)  the  spe¬ 
cific  energy  on  the  initial  isentrope,  and  eH  the  specific  en¬ 
ergy  on  the  Hugoniot.  We  calculated  cv=  1077.0  J  kg-1  K-1 
and  T  =  1 .103  using  published  room-temperature  values  for 
the  heat  capacity  at  constant  pressure,  cp ,30  the  coefficient  of 
volumetric  thermal  expansion,38  isentropic  bulk  modulus, 
Ks,2&  and  ideal  equilibrium  density,  p0. 12  These  values  yield 
AT S=1  A  K.  The  pressure  on  the  isentrope  expanded  to  lead¬ 
ing  order  about  the  equilibrium  specific  volume,  Vo,  is 


Ps(V)  =  Ks 


(4) 


where  Q  is  the  fundamental  derivative.39  The  energy  on  the 
initial  isentrope  can  be  expressed  analytically  as 
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es(V)  =  - 


Ps(V)dV, 


(5) 


An  experimental  value  for  fundamental  derivative  was  ob¬ 
tained  from  the  gradient  of  the  Hugoniot  in  the  Us—  U p  plane 
determined  from  isentropic  compression  experiments  on 
RDX  (100),  namely,  Q= 5.6.39,40  The  energy  on  the  Hugoniot 
is  given  simply  by  the  Hugoniot  equation,21,39 
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Combining  the  values  of  PRH  and  V/V0  from  our  NEMD 
simulations  and  the  calculated  values  of  es(V),  eH,  and  cv 
with  Eq.  (3),  we  estimate  the  additional  heating  arising  only 
from  the  passage  of  the  shock  front,  A7W=289  K.  Hence, 
the  total  temperature  increase  in  real  RDX  during  shock 
loading  is  about  A'/ =363  K.  Thus,  our  classical  MD  simu¬ 
lations  underestimate  the  true  shock  heating  by  factor  of 
about  4,  although  this  should  be  considered  to  be  an  upper 
limit  owing  to  the  temperature  dependence  of  the  experimen¬ 
tal  heat  capacities  and  T. 


3.  Rotational  order 

A  rotational  order  parameter,  P2,  was  used  to  determine 
the  onset  of  shear  band  formation  and  characterize  their  in¬ 
ternal  structure.  P2  characterizes  the  relative  orientation  be¬ 
tween  a  given  vector  within  each  molecule,  in  this  case  the 
vector  connecting  carbon  and  nitrogen  atoms  on  opposite 
sides  of  the  six-membered  ring,  at  time  t=  0  and  that  same 
vector  at  subsequent  time.  P2  is  defined  as 

1  'V  1 

^2(0  =  —2  -[3(ri,.(0)  •  ri,(0)2  -  1],  (8) 

A  ,=1  2 

where  N  is  the  number  of  molecules  and  ri,(f)  the  unit  vector 
between  specified  atoms  in  molecule  i  at  time  t  and  P2  equals 
unity  in  a  system  with  perfect  rotational  order  and  1/4  when 
there  is  no  rotational  order. 

Plots  of  P2  for  a  crystalline  region  and  an  adjacent  shear 
band  in  the  long  and  short  cells  shock  loaded  at  Up 
=  1.0  km  s_l  are  shown  in  Figs.  10(a)  and  10(b),  respec¬ 
tively.  The  rotational  order  parameters  calculated  during  the 
NEMD  and  SFABC  simulations  have  been  merged  so  that 
they  are  plotted  as  a  function  of  the  total  simulation  time.  All 
plots  are  continuous  in  value  and  gradient  across  the  vertical 
lines  which  denote  the  time  at  which  the  SFABC  was  ap¬ 
plied,  providing  convincing  evidence  that  applying  SFABCs, 
in  this  case,  does  not  introduce  unphysical  artifacts  in  the 
simulations. 

Owing  to  thermal  motion,  P2  ~  0.96  in  the  perfect  crystal 
at  300  K.  We  do  not  observe  a  significant  change  in  the  value 
of  P2  after  the  passage  of  a  shock  wave  although  Fig.  7 


FIG.  10.  Rotational  order  parameter  in  a  crystalline  region  and 
shear  band  after  shock  loading  at  Up=  1.0  km  s_1.  The  vertical  line 
denotes  the  time  at  which  SFABCs  were  applied,  (a)  Long  cell  and 
(b)  short  cell. 

shows  notable  differences  in  the  RDFs  of  the  shocked  and 
unshocked  systems.  Hence,  shock  loading  changes  signifi¬ 
cantly  the  crystal  structure  on  a  center-of-mass  level  but  in¬ 
dividual  molecules  largely  retain  their  original  orientation. 

The  rotational  order  parameter  reveals  clearly  the  onset  of 
shear  band  formation  and  its  subsequent  structural  evolution. 
In  both  the  long  and  short  cells,  P2  decreases  rapidly  from 
about  0.96  to  about  0.35  during  the  initial  stage  of  growth 
where  the  rate  of  viscous  heating  is  high  and  amorphization 
progresses  rapidly.  The  initial  stage  of  growth  takes  about 
20-30  ps.  The  shear  bands  subsequently  evolve  relatively 
slowly  to  a  state  of  rotational  disorder  over  a  time  period  of 
about  100  ps.  This  secondary  stage  of  growth  is  associated 
with  the  equilibration  of  the  thermal  gradients  generated  dur¬ 
ing  the  initial  growth  stage,  as  described  in  Sec.  IV  B  2.  The 
rotational  order  parameter  for  the  crystalline  regions  de¬ 
creases  slowly  during  shear  band  formation  from  about  0.96 
to  about  0.88.  We  attribute  this  decrease  to  both  the  40-50  K 
rise  in  temperature  in  the  crystalline  regions  that  results  from 
the  intense  heating  in  the  shear  bands  and  the  misorientation 
of  the  crystalline  regions  with  respect  to  their  original  axes 
due  to  the  shear  strains  mediated  by  viscous  flow  (see  Fig. 

4). 

We  computed  the  autocorrelation  function,  averaged  over 
all  time  origins,  of  the  rotational  order  parameter  for  eight 
representative  molecules  randomly  dispersed  within  a  shear 
band.  By  fitting  the  initial  decay  of  the  autocorrelation  func¬ 
tion  to  C(f)=A  exp (-//  r),  where  A  is  a  constant,  we  estimate 
that  the  molecules  within  the  shear  bands  lose  memory  of 
their  initial  orientation  and  amorphize  within  a  characteristic 
time  interval  of  r*  =  9.7  ps. 

4.  Molecular  conformation 

At  ambient  temperature  and  pressure,  the  point  group  of 
molecules  in  the  a-RDX  crystal  structure  is  Cs.  The  six- 
membered  ring  is  in  a  chair  conformation  with  two  nitramine 
group  N-N  bonds  oriented  axially  (A)  to  the  ring  and  one 
equatorially  (£);  this  molecular  conformation  is  commonly 
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FIG.  11.  Fraction  of  molecules  in  a  given  ring  conformation 
after  shock  loading  at  l/;)=0.63  km  s-1.  The  vertical  line  denotes 
the  time  at  which  SFABCs  were  applied. 

denoted  as  AAE.  Shock  compression  or  melting/ 
amorphization  may  change  the  ring  conformation  and/or  the 
orientations  of  the  three  nitro  groups  with  respect  to  the  ring. 
To  characterize  the  geometry  of  the  molecular  ring  we  cal¬ 
culated  ring-puckering  coordinates  using  the  formalism  of 
Cremer  and  Pople.41  The  orientation  of  the  three  nitro  groups 
with  respect  to  the  ring  was  determined  by  calculating  the 
angles  between  the  three  N-N  bonds  and  the  vector  normal  to 
the  ring. 

In  Fig.  11  we  plot  the  fraction  of  molecules  in  a  given  ring 
conformation  as  a  function  of  the  total  NEMD  plus  SFABC 
simulation  time  during  shock  loading  at  t/p=0.63  kms"1 
(Prh=5.7  GPa).  Before  the  arrival  of  the  shock  wave  in  the 
examined  region  of  the  simulation  cell,  96%  of  the  molecular 
rings  are  in  the  chair  conformation  and  2%  are  in  each  of  the 
half-boat  and  half-chair  conformations.  The  passage  of  the 
shock  wave  through  the  sampled  volume  induces  a  change  in 
the  relative  populations  of  the  three  conformations,  namely, 
the  fraction  of  molecular  rings  identified  as  being  in  the  chair 
conformation  falls  to  85%,  while  those  in  the  half-boat  and 
half-chair  conformations  increase  to  8%  and  7%  of  the  total, 
respectively.  The  fraction  of  molecular  rings  in  the  boat  and 
twist-boat  conformations  is  negligible  both  before  and  after 
the  passage  of  the  shock  front.  Furthermore,  we  did  not  de¬ 
tect  nitro-group  conformations  other  than  AAE.  Since  no  de¬ 
formation  mechanisms  are  active  to  relieve  the  uniaxial 
shock  compression  at  this  Up,  the  relative  populations  of 
molecules  in  these  three  ring  conformations  remain  constant. 

The  same  analysis  was  performed  for  a  crystalline  region 
and  adjacent  shear  band  in  the  long  and  short  cells  shock 
loaded  at  Up=  1.0  km  s"1.  The  fractions  of  molecules  in  a 
given  ring  conformation  as  a  function  of  total  simulation 
time  in  crystalline  regions  in  the  long  and  short  cells  are 
plotted  in  Figs.  12(a)  and  12(b),  respectively.  Behavior  simi¬ 
lar  to  that  observed  under  shock  loading  at  Up 
=  0.63  kms'1  was  found  but  with  some  important  differ¬ 
ences.  First,  the  passage  of  the  shock  front  through  the 
sampled  regions  mediates  a  greater  shift  in  the  relative  popu¬ 
lations  of  the  ring  conformations.  In  both  the  long  and  short 
simulation  cells,  immediately  after  the  passage  of  the  shock 


FIG.  12.  Fraction  of  molecules  in  a  given  ring  conformation  in  a 
crystalline  region  after  shock  loading  at  Up=  1.0  km  s-1.  The  ver¬ 
tical  lines  denote  the  time  at  which  SFABCs  were  applied,  (a)  Long 
cell  and  (b)  short  cell. 

front,  only  60%  of  the  molecular  rings  are  identified  as  chair 
conformers.  The  fractions  of  half-boat  and  half-chair  confor¬ 
mations  do  not  increase  in  equal  amounts;  initially  25%  of 
molecules  adopt  the  half-chair  conformation  and  15%  the 
half-boat  conformation.  As  before,  in  the  crystalline  regions 
we  observe  negligible  fractions  of  the  boat  or  twist-boat  con¬ 
formations  or  nitro  groups  in  orientations  other  than  AAE. 
However,  since  in  these  simulations  at  Up=  1.0  km  s_1  shear 
bands  form  that  relax  the  large  shock-induced  shear  stresses, 
the  crystalline  regions  are  able  to  relax  to  a  state  in  which  the 
average  molecular  conformation  is  more  similar  to  that  in  a 
perfect  a-RDX  crystal.  Hence,  as  shear  bands  grow  through¬ 
out  the  system,  the  fraction  of  molecular  rings  identified  as 
being  in  the  chair  conformation  in  the  crystalline  regions 
increases  to  around  90%  over  100  ps,  while  the  fraction  of 
molecular  rings  in  the  half-boat  and  half-chair  conformations 
decay  accordingly.  Thus,  we  propose  that  plastic  deformation 
via  the  formation  and  growth  of  shear  bands  under  shock 
loading  in  (100) -oriented  RDX  single  crystals  is  largely 
driven  by  the  free  energy  released  upon  relaxation  of  the 
shock-induced  ring  conformations  to  relative  populations 
that  are  more  consistent  with  the  a-RDX  crystal  structure. 

Analyses  of  molecular  ring  conformations  in  the  shear 
bands  in  the  long  and  short  cells  are  presented  in  Figs.  13(a) 
and  13(b),  respectively.  As  for  the  crystalline  regions,  imme¬ 
diately  following  the  passage  of  the  shock  front  around  60% 
of  molecular  rings  are  in  the  chair  conformation  with  25% 
and  15%  in  the  half-chair  and  half-boat  conformations,  re¬ 
spectively.  These  populations  begin  to  relax  initially  due  to 
the  growth  of  shear  bands  elsewhere  in  the  cells.  However, 
once  shear  bands  start  to  propagate  into  these  volumes  and 
melting  commences,  the  populations  of  the  five  possible  ring 
conformations  evolve  to  distributions  more  consistent  with 
the  liquid  state.  In  both  the  long  and  short  cells,  the  fractions 
of  molecular  rings  in  the  chair,  half-boat,  and  half-chair  con¬ 
formations  decrease  over  100  ps  to  34%,  13%,  and  12%, 
respectively.  In  contrast  to  the  crystalline  regions  discussed 
in  the  preceding  paragraph,  20%  and  21%  of  the  molecular 
rings  in  the  shear  bands  are  in  the  boat  and  twist-boat  con- 
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FIG.  13.  (Color  online)  Fraction  of  molecules  in  a  given  ring 
conformation  in  a  shear  band  after  shock  loading  at  Up 
=  1.0  km  s-1.  The  vertical  lines  denote  the  time  at  which  SFABCs 
were  applied,  (a)  Long  cell  and  (b)  short  cell. 

formations,  respectively.  Hence,  the  appearance  of  molecules 
in  these  conformations  signals  the  onset  of  melting  since 
they  serve  as  an  indicator  of  the  increase  in  conformational 
disorder  facilitated  by  the  liquid  state.  Furthermore,  at  the 
same  simulation  times  that  increase  in  the  populations  of 
boat  and  twist-boat  conformations  are  observed,  we  also  find 
notable  increases  in  the  fraction  of  molecules  identified  as 
not  being  in  the  AAE  conformation.  The  fractions  of  mol¬ 
ecules  within  the  shear  bands  in  the  long  and  short  simula¬ 
tion  cells  identified  as  not  possessing  the  AAE 
conformation — that  is,  those  molecules  whose  nitro  groups 
are  in  an  A  EE  or  EEE  conformation — are  plotted  as  a  func¬ 
tion  of  total  simulation  time  in  Figs.  14(a)  and  14(b),  respec¬ 
tively.  In  both  the  long  and  short  cells,  around  30%  of  the 
molecules  in  the  shear  bands  are  not  in  the  AAE  conforma¬ 
tion  once  a  steady-fluctuating  state  is  achieved.  Although  the 
AAA  conformation  is  suggested  to  be  favored  in  the  gas 
phase,42'43  we  detect  only  a  negligible  number  of  molecules 
in  this  conformation  in  the  condensed  phase. 


FIG.  14.  Fraction  of  molecules  in  a  shear  band  with  nitro  groups 
not  in  the  AAE  conformation  after  shock  loading  at  Up 
=  1.0  km  s-1.  The  vertical  lines  denote  the  time  at  which  SFABCs 
were  applied,  (a)  Long  cell  and  (b)  short  cell. 
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V.  DISCUSSION  AND  CONCLUSIONS 
A.  Efficacy  of  the  SFABC  method 

The  use  of  SFABCs  allowed  us  to  extend  simulation  times 
to  the  extent  that  shock-induced  shear  bands,  which  evolve 
slowly  on  the  time  scale  of  a  typical  NEMD  simulation, 
could  be  studied  over  their  complete  evolution  to  a  steady- 
fluctuating  state.  This  is  demonstrated  with  greatest  clarity  in 
Figs.  13  and  14  where  relative  abundances  of  the  ring  and 
nitro-group  conformations  at  the  point  of  maximum  com¬ 
pression  in  the  respective  NEMD  simulations,  as  denoted  by 
the  vertical  broken  lines,  differ  significantly  from  those  once 
a  steady  state  is  reached.  More  generally,  the  absence  in  Figs. 
10-14  of  perturbations  associated  with  the  application  of  the 
SFABC  at  the  moment  of  maximum  compression  further 
supports  the  validity  of  the  approach. 

There  are  certain  situations  where  the  application  of 
SFABCs  may  not  be  appropriate.  For  example,  it  would  be 
difficult  to  apply  this  method  in  the  study  of  shock  waves  in 
ductile  metals  where,  in  the  two-wave  regime,  an  elastic  pre¬ 
cursor  may  lead  a  strong  plastic  wave  by  tens  of  nanometers. 
In  this  case,  if  SFABCs  are  applied  when  the  elastic  precur¬ 
sor  reaches  the  free  surface,  the  slower  plastic  wave  may 
reflect  off  the  second  piston,  setting  up  undesirable  density 
waves  in  the  material.  However,  this  problem  could  be  over¬ 
come  within  the  SFABC  framework  by  applying  the  second 
piston  within  the  simulation  cell  at  a  point  that  has  been 
passed  by  both  the  elastic  and  plastic  waves  rather  than  at  the 
end  of  simulation  cell.  This  approach  would  still  require  a 
computationally  expensive  NEMD  simulation,  but  subse¬ 
quent  MD  simulations  could  be  performed  using  a  smaller 
number  of  atoms.  In  the  case  of  shock  waves  in  RDX,  we  do 
not  find  strong  plastic  waves  that  lag  significantly  behind  the 
elastic  precursor;  hence  the  SFABC  method  is  particularly 
adept  at  capturing  the  simulation  cell  at  maximum  compres¬ 
sion. 

Hugoniostat  methods  have  been  developed  that  allow 
shock  loading  to  be  simulated  within  the  framework  of  MD 
and  which  are  less  computationally  demanding  than  large- 
scale  NEMD  simulations.19'44  47  These  methods  modify  the 
equation  of  motions  of  atoms  to  drive  a  system  toward  a 
prescribed  state  on  the  Hugoniot  locus.  It  is  not  obvious  that 
Hugoniostats  are  appropriate  tools  for  use  in  the  simulation 
of  molecular  crystals  since  phenomena  at  the  shock  front  in 
these  materials  are  extremely  complex  owing  to  their  rela¬ 
tively  low  symmetry  and  many  intramolecular  degrees  of 
freedom.  Explicit  simulation  of  the  shock  front  using  NEMD 
will  provide  a  more  accurate  description  of,  for  example,  the 
initial  overshoot  in  intermolecular  temperature  and  its  equili¬ 
bration  with  the  intramolecular  vibrational  modes  in  the  na¬ 
nometers  behind  the  shock  front. 14  27  Previous  work  showed 
that  these  phenomena  are  of  fundamental  importance  to  the 
homogeneous  nucleation  of  dislocations  in  (lll)-oriented 
RDX  single  crystals.14  However,  Hugoniostats  do  provide 
certain  advantages  over  SFABCs;  for  example,  the  former 
ensures  that  the  system  remains  on  the  Hugoniot  locus,  while 
the  latter,  in  general,  does  not.  SFABCs  capture  the  simula¬ 
tion  cell  at  a  constant  volume  that  corresponds  only  to  the 
system  after  the  passage  of  the  leading  shock  wave(s). 
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Shock-induced  phase  transformations  or  plastic  deformation 
that  lags  the  elastic  wave  will  reduce  the  diagonal  component 
of  the  stress  tensor  parallel  to  the  direction  of  shock  propa¬ 
gation  and  the  system  under  SFABCs  will  depart  from  the 
Hugoniot  locus.  Hugoniostat  methods  enforce  a  constant 
pressure  rather  than  a  constant  volume  simulation.  This  limi¬ 
tation  of  the  SFABC  method  may  be  acceptable  in  certain 
circumstances,  for  example,  if  one  wishes  to  extend  simula¬ 
tion  times  to  identify  events  that  occur  at  a  rate  too  low  to  be 
observed  directly  in  NEMD  simulations,  such  as  the  homo¬ 
geneous  nucleation  of  dislocations,  second  phases,  or  the  ini¬ 
tiation  of  slow,  complex  chemical  reactions.  Hence,  the 
SFABC  method  may  be  of  particular  advantage  in  cases 
where  the  description  of  interatomic  bonding  is  sufficiently 
expensive  that  large-scale  NEMD  simulations  are  not  fea¬ 
sible  but  an  accurate  treatment  of  the  dynamics  associated 
with  the  shock  front  is  required.  The  principal  advantages  of 
SFABCs  over  Hugoniostats  are  that  the  former  are  extremely 
simple  to  implement,  robust,  do  not  affect  the  equations  of 
motion  of  atoms  artificially,  and  provide  an  essentially  exact 
treatment  of  phenomena  at  the  shock  front. 

In  order  to  reach  the  same  total  NEMD  plus  SFABC 
simulation  time  for  the  material  subvolumes  sampled  here  in 
a  normal  NEMD  simulation,  a  simulation  cell  of  length 
604a,  containing  over  40X106  atoms,  would  be  required. 
Such  simulations  would  be  computationally  challenging  at 
present  (the  simulation  reported  here  using  the  quasi-2D  cell 
required  8.8  CPU  years).  Thus,  the  SFABC  method  facili¬ 
tates  simulations  that  otherwise  may  be  prohibitively  expen¬ 
sive.  Furthermore,  the  results  reported  in  Sec.  IV  B  for  the 
long  and  short  cells,  the  latter  containing  half  the  number  of 
atoms  of  the  former,  concur  in  every  regard  once  a  steady- 
fluctuating  state  is  achieved  (these  simulations  appear  to  dif¬ 
fer  only  in  the  rate  of  shear  band  formation).  Thus,  SFABCs 
not  only  enable  simulation  times  to  be  extended  for  a  given 
cell  size  but  also  facilitate  simulations  with  significantly 
fewer  atoms  with  no  loss  in  the  accuracy  of  the  description 
of  the  final  state.  Finally,  in  comparison  with  Hugoniostat 
methods,  use  of  the  SFABC  provides  an  exact  treatment  of 
the  dynamics  during  shock  loading  and  is  far  more  simple 
and  robust  in  its  implementation. 

B.  Implications  for  initiation  sensitivity 

Dislocation-mediated  plasticity  was  not  observed  during 
shock  loading  at  either  f/p=0.63  km  s_1  (PRH=5.7  GPa)  or 
Up=  1.0  kms-1  (Prh=9.7  GPa).  This  result  was  expected 
for  this  orientation  of  the  shock  propagation  direction  since 
none  of  the  slip  systems  identified  under  quasistatic  loading 
are  subject  to  a  resolved  shear  stress.1’  Furthermore,  we  can 
eliminate  the  possibility  that  for  this  orientation  of  the  shock 
propagation  direction,  other  previously  unknown  slip  sys¬ 
tems  are  activated  by  the  passage  of  the  shock  wave;  that  is, 
phenomena  such  as  those  observed  during  shock  loading  on 
(111)  (Ref.  14)  do  not  occur  for  this  orientation.  Thus,  based 
on  the  steric -hindrance  model,  we  expect  high  sensitivity  for 
impacts  on  (100),  in  accord  with  Ref.  13. 

An  unanticipated  deformation  mechanism  was  identified 
during  shock  loading  on  (100)  at  Up=l.O  kms'1,  namely. 


the  formation  and  growth  of  shear  bands.  While  the  steric- 
hindrance  model  suggests  that  deformation  mechanisms 
which  serve  to  relax  shear  stresses  promote  low  impact  sen¬ 
sitivity,  the  intense  heating  caused  by  viscous  flow  within  the 
shear  bands  will  certainly  promote  thermal  molecular  de¬ 
composition.  Moreover,  viscous  flow  may  promote  the  me¬ 
chanical  rupture  of  intramolecular  bonds,  further  enhancing 
the  rate  of  molecular  decomposition.  Thus,  in  a  sense,  shear 
bands  can  be  considered  to  be  homogeneously  nucleated  hot 
spots  since  the  potential  energy  stored  in  the  shock  com¬ 
pressed  material  is  partly  released  in  spatially  localized  ma¬ 
terial  subvolumes. 

One  might  expect  that  if  dislocation  slip  systems  could  be 
activated  at  relatively  low  values  of  PRH  at  this  orientation  of 
the  shock  propagation  direction,  then  shear  stresses  would 
not  build  to  the  level  where  deformation  via  shear  bands  is 
inevitable.  However,  NEMD  studies  of  shock  waves  in 
ff-HMX  by  Jaramillo  et  al r1  revealed  a  gradual  transition 
from  dislocation-mediated  plasticity  to  deformation  mediated 
by  the  growth  of  liquidlike  shear  bands  with  increasing  PRH. 
Hence,  the  deformation  of  energetic  molecular  crystals  via 
the  nucleation  and  growth  of  shear  bands  appears  to  be  re¬ 
lated  to  intrinsic  instabilities  of  the  crystal  structure. 

C.  Comments  on  the  effects  of  nonclassical  phenomena  in 
simulations  of  nonequilibrium  processes  in  molecular 
crystals 

It  is  well  known  that  in  the  temperature  range  of  interest 
in  this  study,  the  heat  capacity  of  a  system  of  classical  oscil¬ 
lators  with  a  frequency  distribution  representative  of  a  mo¬ 
lecular  crystal  is  significantly  larger  than  the  value  in  the 
corresponding  quantum-mechanical  system;  thus,  for  the 
weak-to-moderate  shock  strengths  considered  here,  the  shock 
temperature  predicted  on  the  basis  of  classical  MD  will  be 
lower  than  the  “real”  value,  as  was  shown  in  Sec.  IV  B  2. 
Nevertheless,  this  underestimate,  which  is  inherent  to  the  use 
of  classical  mechanics,  does  not  affect  our  underlying  theses, 
namely,  that  the  homogeneous  nucleation  of  liquidlike  shear 
bands  occurs  in  [100]-oriented  RDX  crystals  under  suffi¬ 
ciently  strong  shock  loading  and  the  viscous  flow  of  material 
within  these  defects  leads  to  an  additional,  non-negligible 
localized  heating.  Since  our  simulations  underestimate  the 
magnitude  of  the  shock  heating,  it  is  a  reasonable  assumption 
that  experimentally  these  defects  will  nucleate  at  signifi¬ 
cantly  lower  shock  pressures.  Furthermore,  the  magnitude  of 
the  heating  in  the  shear  bands  will  be  higher  than  the  70-90 
K  we  calculate  since  the  “error”  in  the  heat  capacity  of  ma¬ 
terial  due  to  classical  mechanics  is  equally  present  within 
these  subvolumes.  Thus,  while  our  calculated  shock  and  de¬ 
fect  temperatures  are  low  compared  to  those  in  real,  quantum 
dynamical  RDX,  they  are  almost  directly  proportional;  that 
is,  based  on  specific  heat  alone,  classical  MD  simulations 
will  also  underestimate  the  additional  heating  within  the 
shear  bands  by  roughly  the  same  factor  of  4  as  was  predicted 
for  the  shock  heating. 

Fundamentally,  the  “specific-heat  problem”  of  interest 
here  for  the  case  of  a  condensed  phase  molecular  system  and 
the  “zero-point  energy  problem”  that  arises  for  the  case  of 
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gas-phase  polyatomic  molecules  are  manifestations  of  quan¬ 
tum  effects  and  have  no  classical  analog.  At  present  there  is 
no  satisfactory  method  for  performing  large-scale  MD  simu¬ 
lations  with  a  quantum-mechanically  accurate  treatment  of 
the  zero-point  energy  problem  (for  gas-phase  species)  or  the 
specific-heat  problem  (for  condensed  molecular  phases). 
Nevertheless,  our  simulations  were  designed  to  mitigate  in  a 
physically  sensible  way  the  overestimation  error  in  the  spe¬ 
cific  heat  of  RDX  by  freezing  out  via  geometric  constraints, 
the  contributions  to  the  heat  capacity  from  vibronic  degrees 
of  freedom  with  the  highest  modal  Debye  temperature, 
namely,  the  C-H  bond  stretches.  Attempts  to  go  beyond  this 
type  of  strategy  are  problematic:  Thompson  and 
co-workers48,49  demonstrated  that  for  model  Hamiltonians 
with  small  numbers  of  degrees  of  freedom  that  both  “active” 
and  “passive,”  quantum-inspired  corrections  to  classical  dy¬ 
namic  trajectories  may  result  in  severely  unphysical  behav¬ 
ior.  Specifically,  even  for  systems  containing  only  zero-point 
vibrational  energy,  observables  obtained  from  uncorrected 
classical  MD  trajectories  resembled  more  closely  the  true 
quantum-mechanical  ones  than  those  obtained  for  trajecto¬ 
ries  that  included  quantum-corrected  classical  dynamics. 

Another  option,  essentially  a  generalization  of  the  idea  of 
freezing  out  the  high-frequency  degrees  of  freedom  em¬ 
ployed  in  the  current  study,  is  to  introduce  a  Debye-type 
dependence  of  the  heat  capacity  on  temperature  using  the 
mesoparticle  dynamics  of  Strachan  and  Holian.50  In  this  ap¬ 
proach,  a  reduced-dimension  representation  of  the  potential- 
energy  function  results  in  a  significant  decrease  in  the  clas¬ 
sical  heat  capacity  of  the  system.  This  is  balanced  by 
introducing  into  the  equations  of  motion  for  the  system  “im¬ 
plicit  degrees  of  freedom”  that  correspond  to  local  heat  res¬ 
ervoirs  within  the  mesoparticles.  These  internal  degrees  of 
freedom  can  be  populated  using  either  classical  or  quantum 
statistical  descriptions  of  their  behavior.  While  this  approach 
is  attractive,  it  is  rather  empirical  and  mesopotentials  pres¬ 
ently  do  not  demonstrate  the  levels  of  sophistication  needed 
to  describe  accurately  a  system  as  complex  as  an  RDX  crys¬ 
tal  under  extreme  loading  conditions.  Thus,  although  simu¬ 
lations  based  on  classical  MD  display  certain  well-known 
and  well-understood  limitations,  at  present  there  exists  no 
better  tool  for  the  study  of  shock-induced  phenomena  in  or¬ 
ganic  molecular  crystals. 
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APPENDIX:  CONSERVATION  OF  ENERGY  IN  SFABC 
SIMULATIONS 

As  discussed  briefly  in  Sec.  I,  shock  waves  can  be  gener¬ 
ated  in  NEMD  simulations  by  driving  a  piston  at  a  specified 
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Up  into  a  stationary  simulation  cell  as  in  Ref.  22  or  by  im¬ 
pacting  a  simulation  cell  onto  a  fixed  rigid  piston  at  a  speci¬ 
fied  Up.  For  convenience  and  clarity,  we  refer  hereafter  to 
these  approaches  as  the  moving-piston  and  moving-cell 
methods,  respectively.  These  methods  are  Galilean  invariant, 
thus  it  is  largely  a  matter  of  personal  preference  which  is 
employed.  However,  we  found  unexpected  differences  be¬ 
tween  these  formally  equivalent  methods  upon  the  applica¬ 
tion  of  SFABCs. 

If  the  moving-piston  method  is  used  to  generate  the  shock 
wave  in  a  NEMD  simulation,  SFABCs  can  be  applied  in 
much  the  same  way  as  described  in  Sec.  III.  The  time  at 
which  the  simulation  cell  reaches  maximum  compression  can 
be  determined  with  high  accuracy  by  again  plotting  as  a 
function  of  time  the  x  coordinate  of  the  center  of  mass  of  the 
last  molecule  and  identifying  when  the  curvature  becomes 
nonzero.  The  second  piston  is  then  applied  by  setting  the 
force  and  velocity  vectors  of  all  atoms  within  a  specified 
distance  from  the  end  cell  to  zero  and  (Up, 0,0),  respectively. 
In  this  manner,  the  volume  of  the  cell  is  fixed  to  that  corre¬ 
sponding  to  maximum  compression  and  the  entire  system 
translates  with  a  velocity  (Up, 0,0). 

When  a  shock  wave  was  generated  in  (100) -oriented  RDX 
using  the  moving-cell  method  at  Up=  1.0  kms'1  and 
SFABCs  were  applied  as  in  Sec.  Ill,  during  equilibrium  MD 
simulations  in  the  microcanonical  ensemble,  total  energy  was 
conserved  to  a  tolerance  of  1.2  ppm.  However,  when  the 
moving-piston  method  was  employed  we  found  that  the  total 
energy  was  conserved  to  only  2.2  parts  per  thousand  using 
otherwise  identical  computational  protocols.  Despite  the 
thousandfold  decrease  in  energy  conservation  in  the  latter 
simulation,  we  found  that  the  relative  particle  trajectories  in 
the  two  simulations  were  identical  until  numerical  diver¬ 
gences  arose  due  to  chaotic  dynamics. 

The  origin  of  the  apparent  poor  conservation  of  the  total 
energy  when  the  simulation  cell  translates  is  determined  en¬ 
tirely  by  the  contribution  from  the  kinetic  energy,  T,  to  the 
total  energy.  We  can  write  the  kinetic  energy  of  an  ensemble 
of  N  thermalized  particles  of  mass  m;  in  a  system  whose 
center  of  mass  translates  with  constant  velocity  U 
=  (Up, 0,0)  as 

1  N 

7TU)  =  -2  K-V;  •  V;  +  m,U  •  U  +  2m, V;  •  U],  (Al) 

2  ,=i 

where  ,v? ,vf)  is  the  velocity  of  particles  with  respect 

a  coordinate  system  translating  with  velocity  U.  The  first 
term  is  related  to  the  temperature  of  the  ensemble  and  the 
second  term  is  a  constant.  The  third  term,  T ,  is  a  contribu¬ 
tion  to  the  kinetic  energy  arising  from  the  momentum  of  the 
system  parallel  to  U, 

N 

T  =  UrX  mtf.  (A2) 

1=1 

Hence,  the  apparent  poor  conservation  of  the  total  energy  of 
simulation  cells  under  SFABCs  translating  with  velocity  U p 
is  explained  by  fluctuations  of  the  total  momentum,  1,f=1  mp*. 
Hence,  in  the  microcanonical  ensemble,  MD  configuration 
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FIG.  15.  (a)  Oscillation  of  the  total  momentum  in  a  system 
under  SFABCs  translating  with  velocity  at  Up=  1.0  km  s-1  and  (b) 
deviations  of  the  total,  kinetic,  and  potential  energies  of  the  system 
from  their  mean  values. 


trajectories  are  independent  of  U  (i.e.,  they  are  Galilean  in¬ 
variant  in  configuration  space),  but  the  kinetic  energy  of  the 
system  depends  on  the  scalar  product  of  the  total  momentum 
with  U. 

To  illustrate  the  contribution  from  the  nonconserved  linear 
momentum  to  the  kinetic  energy  in  a  translating  system  we 
plot  in  Fig.  15(a)  the  quantity  m-v)  as  a  function  of  simu¬ 
lation  time  for  an  RDX  crystal  under  SFABCs.  The  simula¬ 
tion  cell  was  shock  loaded  using  the  moving-piston  method 
and  is  translating  with  velocity  Up=  1.0  km  s-1.  In  Fig.  15(b) 
we  plot  as  a  function  of  simulation  time  the  deviation  of  the 
total,  kinetic,  and  potential  energies  of  the  system  from  their 
mean  values  time  averaged  over  the  simulation.  Comparing 
Figs.  15(a)  and  15(b)  it  is  evident  that  the  kinetic  and  total 
energies  track  oscillations  of  the  total  momentum,  while  the 
potential  energy,  which  reflects  the  relative  particle  trajecto¬ 
ries  of  the  system,  is  independent  of  the  total  momentum. 
The  slow  decrease  in  the  potential  energy  in  Fig.  15(b)  arises 
from  the  nucleation  and  growth  of  shear  bands  in  the  system. 
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